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ABSTRACT 

A fully Sinc-Galerkin method for recovering the spatially varying stiffness parameter 
in fourth-order time-dependent problems with fixed and cantilever boundary conditions is 
presented. The forward problems are discretized with a sine basis in both the spatial and 
temporal domains. This yields an approximate solution which converges exponentially and 
is valid on the infinite time interval. When the forward methods are applied to parameter 
recovery problems, the resulting inverse problems are ill-posed. Tikhonov regularization is 
applied and the resulting minimization problems are solved via a quasi-Newton/trust region 
algorithm. The L-curve method is used to determine an appropriate value of the regulariza- 
tion parameter. Numerical results which highlight the method are given for problems with 
both fixed and cantilever boundary conditions. 


1 This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 


l 



1 Introduction 


In this paper, a fully Sinc-Galerkin method is introduced for the numerical recovery of mate- 
rial parameters in fourth-order time-dependent problems. To illustrate the method, consider 
the problem of estimating the spatially varying parameter EI(x ) in the state equations 

«">* -Sr+aK^S?)"**’ 0, 0<I<I ‘ >0 

u(0, t) = ti(l, t) = 0, t> 0 (1.1) 

d £(o, t) -£(M)-o, <>o 

du 

u(x,0)= — (x,0) = 0, 0 < x < 1 

and 

C(EI)u = /(x,<), 0 < x < 1, t > 0 

u(0,4)=5((), = m t>0 (1.2) 

a £(o,t)=m, ‘>° 

u(x, 0) = — (x,0) = 0, 0 < x < 1 

l/t 

given measurements of the data at the points {(x p , ^ 9 )}p=l|...’n* i n (0, 1) x 2R + . These formula- 
tions are generalizations of the equations which arise when using the Euler-Bernoulli theory 
to model beams with flexural rigidity EI(x) and fixed and cantilever ends, respectively. For 
ease of presentation throughout the paper, the boundary conditions in (1.1) and (1.2) will 
be referred to as fixed and cantilever conditions with the general y(t) and ?(t) included to 
allow for boundary controllers. 

Since EI(x ) denotes the flexural stiffness, it is physically reasonable to assume that El 
is continuous on [0, 1] and to let the admissible parameter Bet Q be defined by 

Q = {El € H\ 0, 1) : EI(x) > EI 0 > 0} 

(see (9J). With this definition, the existence of a unique solution u to the forward problem 
can be obtained on any fixed time interval [0,r],r > 0, for / sufficiently smooth. 
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In order to apply the results from classical operator theory, the inverse problems corre- 
sponding to (1.1) and (1.2) can be written as operator equations of the form fC(EI) = d 
where K is compact. The procedure for problems with fixed boundary conditions is out- 
lined below with the formulation for problems with cantilever boundary conditions following 
similarly. Further theory for this latter case can be found in [1]. 

In formulating the operator equation corresponding to (1.1), it is easiest to first consider 
the spatial problem 


( EI(x)u — q(x)u = /(*), 0 < x < 1 

u(0) = u(l) = tt'(O) = u'(l) = 0 


(1.3) 


where EI(x ) and q(x) are both strictly positive on [0, 1]. From the theory of [10] as noted 
in [20], there exists a Green’s function for (1.3) which will be denoted by G[x, s, EI(s)) to 
emphasize its dependence on El. It follows that G and are continuous on [0,1] and that 
the state solution u is given by 


u(x) = f G(x, s, EI(s))f(s)ds. 
Jo 


(see [20], pages 205-206). 

These one-dimensional results can then be extended to the time-dependent problem (1.1) 
via a separation of variables. The substitution of u(x, t) — X(x)T(t) into the homogeneous 
problem corresponding to (1.1) yields 


T'\t) + \T(t) = 0 


and the eigenvalue problem 


(EI(x)X"(x))" - \X(x) = 0, 0 < x < 1 

X(0) = X(l) = X'(0) = X't 1 ) = o 

where A > 0. From arguments similar to those in [5] and [20], it follows that since EI{x) > 0 
on [0,1], the state solution to (1.1) can be represented by 

u(x, t) = f Q(x, t, s, EI(s))d3. 


2 



The function Q depends on a Green’s function as well as the expansions of T(t) and the forcing 
function /. Hence Q is continuous with respect to El and £§j is bounded on (0, 1) x (0, r] 
for / sufficiently smooth. 

The inverse problem can then be written as the operator equation 

K{EI) = d (1-4) 

where d denotes the data. The nonlinear operator JC : H 2 ( 0, 1) — ► I> 2 ((0, 1) x (0, r]) is defined 

by i 

K{EI) = C [' £(•, •, s, EI(s))ds (1-5) 

Jo 

where the observation operator C maps the state solution into the infinite dimensional ob- 
servation space by 

Cip{x,t) = {tp(x,t)} ; (1-6) 

that is, C samples functions continuously throughout (0,1) x (0,r). Note that the choice 
L 2 ((0, 1) x (0, t]) for observation space is physically reasonable. 

To show that K. is compact, it is first shown that it is weakly continuous; that is, 
EI n A El in H 2 ( 0,1) implies that K{EI n ) -+ K.{EI) in the L 2 norm. For every (x, t) 
in (0, 1) x (0, r] it follows that 

I cf G(x,t,s,EI n (s))ds-C [ £/(x, t, s, EI(s))ds 

1 JO Jo 

< *■'■«*) l^n(») - 

<M I' |E/„( S ) - EI(s)\ds 

Jo 

< M max \EI n (a) — jE7(s)| 

«G[0,1) 

with the last inequality resulting from the continuity of El on [0,1]. Since weak conver- 
gence in 77 3 (0, 1) implies uniform convergence, it follows that K{EI n ) converges uniformly 
to K{ El). The weak continuity of K results from the fact that the uniform norm is stronger 
than the L 2 norm. Hence the operator 1C as defined in (1-5) is compact since weak continuity 
implies compactness. 

Although the procedure just outlined was for the problem (1.1), similar results can be 
obtained for (1.2) once a Green’s function has been found which satisfies the boundary 
conditions. Further theory on problems of this type can be found in [1] and [10]. 
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Consider now the well-posedness of (1.4). First, since Q is continuous, the range of 
the operator 1C lies in ^((0,1) x (0,r]) for any El G Q. Hence there exist elements 
d G L 2 ((0, 1) x (0,t]) for which (1.4) has no solution. Furthermore, since fC is a compact 
operator with an infinite dimensional range, it follows that the Moore-Penrose generalized 
inverse is discontinuous. This in turn indicates that small perturbations in the data d 
may give rise to arbitrarily large perturbations in the solution El 6 Q. Consequently, some 
sort of regularization (i.e., stabilization) is required to obtain an accurate approximation for 
El. 

The regularization technique that is used is Tikhonov regularization [24] , and the problem 
(1.4) is replaced by the minimization problem 

minr„(E/) (1.7) 

where 

%{EI) = j{||K(E/) - d||’ + aJ(EJ)}. 

Here a > 0 is a regularization parameter which controls the tradeoff between goodness of fit 
to the data and stability. The penalty functional J(EI) provides stability and allows the 
inclusion of a priori information about the true parameter EL Since El is assumed to be 
“smooth” in the sense that El £ i/ 2 (0, 1), the penalty functional is taken to be the norm 

J{EI) = \\EI\\ 2 q = j\EI"{x)] 2 dx + e [\eI(x)] 2 <Lx. (1.8) 

JO JO 

with e of order 10 -6 . The reasons for including the second term and forcing J to be strictly 
positive will be discussed in the fourth section of the paper. By using arguments similar to 

those in [7] and [16] and assuming that IC(EI) is one to one, it can be shown that with this 

definition for J(E1) } the solutions EI a to (1.7) converge as the regularization parameter 
a-+0 and as the perturbations in the data and operator tend to zero. 

Due to the infinite dimensionality of Q and that of the state space, the problem (1.7) is an 
infinite dimensional minimization problem. In order to develop a practical numerical scheme, 
the problem must be replaced by a sequence of finite dimensional problems; that is, one 
must approximate the operator K. and minimize the functional T a over a finite dimensional 
admissible subspace of Q. 
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The evaluation of K,(EI ) requires the solution of the partial differential equations (PDE’s) 
(1.1) or (1.2). Similar PDE’s must be solved to obtain the components of the derivative 
K\EI). The construction of an approximate solution to these forward problems commonly 
begins with a Galerkin discretization of the spatial variable with time-dependent coefficients. 
This yields a system of ordinary differential equations which is solved via differencing tech- 
niques. Due to stability constraints on the discrete evolution operator, low order methods 
with small time steps are often required to obtain accurate approximations. Moreover, this 
time-stepping must be repeated at each step in the minimization of (1.7). A final difficulty 
lies in the need to interpolate at data points which do not coincide with the nodes of the 
ODE solver. 

In contrast, the method of this work implements a Galerkin scheme in time as well as 
space. This method thus bypasses many of the difficulties associated with time-stepping 
methods in the context of inverse problems. Corresponding results for the heat equation can 
be found in [12] and [19]. 

The fully Sinc-Galerkin method in space and time has many salient features due both 
to the properties of the basis functions and the manner in which the problem is discretized. 
Perhaps the most distinctive feature of the method is the resulting exponential convergence 
rate when solving the corresponding forward problems. Furthermore, the judicious choice 
of a conformal map provides approximate solutions which are valid on the infinite time 
interval rather than only on a truncated time domain. Finally, the discrete system requires 
no numerical integrations to fill either the coefficient matrices or the right-hand side matrix. 
All three features prove to be advantageous when solving the forward problems and hence 
the inverse problem. 

The foundations of the Sinc-Galerkin method are described in Section 2. The fundamental 
quadrature rules are given, and the exponential convergence rate of this method is stated. 
A thorough review of sine function properties can be found in [22] and [23]. 

In the third section, the Sinc-Galerkin systems for the forward problems are constructed 
and implementation details are discussed. The section includes the outline of a very robust 
and accurate algorithm for solving the resulting matrix systems. 

Section 4 includes the finite dimensional minimization problem with the discussion cen- 
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tering around the construction of the various components of the Tikhonov functional. The 
resulting unconstrained optimization problems are solved via a quasi-Newton/trust region 
algorithm as described in [2] and [8]. 

Numerical results are presented in Section 5. Of the many examples tested, those dis- 
cussed in this section best exhibit the features necessary for the practical implementation 
of the method. A brief discussion of the L - curve technique [6] for determining the regular- 
ization parameter a is given at the beginning of the section, and the applicability of this 
technique in conjunction with the Sinc-Galerkin method is demonstrated by the numerical 
results. Finally, results are included both from data sets with white noise and from data sets 
to which no noise was added. As shown in these examples, the Sinc-Galerkin method works 
equally well in both cases. 


2 Sine Function Properties 


For the Sinc-Galerkin method, the basis functions are derived from the Whittaker cardinal 
(sine) function 

sin(7rx) 

, — oo < x < oo 

7T X 


( 2 . 1 ) 


sinc(x) 


and it translates 

-kh\ 

— )■ h>0 ■ 

For h* = three adjacent members of this sine family (S(k,h*)(x),k = —1,0,1) are shown 
in Figure 1. 


S(k , h)(x ) = sine 


( 
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Figure 1. Three Adjacent Members h*)(x), k - -1,0, 1, h* = of the Translated Sine 
Family. 


To construct basis functions on the intervals (0,1) and (0,oo), respectively consider the 
conformal maps 



* W = /n (l-z) 

(2.2) 

and 


T(u>) = /n(u>). 

(2.3) 

The map <f> carries the eye-shaped region 



C E ={z = l + .»:|<.rj( 1 _ J|<rf< 2 } 

(2.4) 

onto the infinite strip 


D s = {£ = C + *V ’ M < d < ^}- 

(2.5) 

Similarly, the map T 

carries the infinite wedge 



Dw = = t + is : |ar^(iw)| < d < 

(2.6) 


onto the strip These regions are depicted in Figure 2. 
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Figure 2. The Domains Ds, De, and Dw- 


The sine gridpoints z k € (0, 1) in De will be denoted x k since they are real. Similarly, the 5 
gridpoints w k G (0,oo) in Dw will be denoted t k . Both are inverse images of the equispaced 
grid in Ds', that is, 

*» = t-'(kh) = I -^SJ (2 7) 

and 


t k = T~ 1 (kh) = e kh . 


( 2 . 8 ) 


To simplify notation throughout the remainder of this section, the pairs <f > , De and T, Dw = 
are referred to generically as x> D. It is understood that the subsequent definition, theorems, 
and identities hold in either setting. Furthermore, the inverse of x is denoted by ip. 

The important class of functions for sine interpolation and quadrature is denoted B(D) 
and defined next. 
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Definition 2.1. Let B{D ) be the class of functions F which are analytic in D, satisfy 

f \F(z)dz\ — ♦ 0, t — ► ±oo 

J<4>(t+L) 

where L = {is : |s| < d < =}, and on the boundary of D ( denoted dD) satisfy 

N(F) = / \F{z)dz\ < oo. 

J QD 

The following theorem for functions in B{D ) is found in [21]. 

Theorem 2.1. Let T be (0,1) or (0,oo) when x = <f> or T, respectively. If F € B(D) and 
Zj = if(jh) = J = 0, ±1,±2, - • • , then for h > 0 sufficiently small 

(2 - 9) 

Theorem 2.1 illustrates the exponential convergence rate which is a trademark of the sine 
methods. There is a common occasion when it is possible to evaluate the infinite series 
appearing in (2.9), namely when integrating against S(k } h) o x • In general, however, the 
series must be truncated. With additional hypotheses, it is proven in [11] and [22] that the 
truncation need not be at the expense of the exponential convergence. 


Theorem 2.2. Assume F G B(D ) and that there exist positive constants K,a, and fi such 


Fir) 

X'( T ) 

Then for h sufficiently small 


e -a| X (r)| } T G V’(( — OOj 0)) 
T £ ^»([0,oo)). 


( 2 . 10 ) 


/ F(z)di - h £ < K,e- M ' h + ~e~ aMh + (211) 

Jr ' ’ yirW x'M a p 


Theorems 2.1 and 2.2 are used to establish a uniform error bound when constructing an 
approximate solution to the forward fourth-order time-dependent problems. The application 
of these quadrature theorems is facilitated by the identities 


= P(p, h ) 0 xWl 

z- 


1, i = p 
0 | * 7 ^ Py 


( 2 . 12 ) 
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«>■* 


d X 


S(p,h)o X (z) 


Z=Zi 


0, i = p 

(— l)* -p 


(*“P) 


» * ^ P> 


(213) 


S<? = fc 1 


<? 


dx : 




3 ’ 


z = p 


(-2K-1)-' 

(.-p)’ ' ’* p ' 


4 3) ■ 


<z 3 


L rf x 


S(p, /i)oxW 


= < 




0, 

(-ir p 

l (*-p) 3 


z = p 

[6 — 7r a (i — p) a ], i^p, 


and 


(2.14) 


(2.15) 


4!> = h< 

P* 


d* 


d X 


- i S(p,h)ox(z) 


7T 

¥’ 


i = p 


~^ _p). [6 -»*(«-?)*], i^P 


(2.16) 


which denote the evaluation at the gridpoint Z{ of the sinc-map compositions and their 
derivatives with respect to the map x ■ 


3 The Forward Problem 


Two forward problems of interest are 

d 7 u d 2 ( d 2 u \ 

Cu(x,t) = —(x,t) + —{EI(x)-^(x,t)\=f(x,t), 0 < i < 1 <>0 

u(0, t) = u(l, t) = 0, t>0 (3.1) 

s(o.O-|(i.«)-o, «>o 

du 

u(x, 0) = — (x,0) = 0, 0 < x < 1 
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and 


Cu(x,t) = f(x,t), 0 < x < 1, t> 0 

u(o ,<) = «(*), (M) = 7(0. *>° ( 3 - 2 ) 

|(0. o-wo. 

u(i,0) = | r (x,0) = 0, 0<x<l. 

Since a thorough derivation of the Sinc-Galerkin method for problems of this type is 
given in [18], the following discussion contains only that material which is needed for the 
construction of the associated matrix systems. 

To define the Sinc-Galerkin approximation to (3.1), let Si(x) = S{i,h x ) o ^(x) and 
s;(t) = S(j, h t ) o T(f), and take the basis to be where 

s h (x,t) a Si(x)s;(t). 

The approximate solution is then defined by way of the tensor product expansion 

N m N t m x = M x 4- N x 4- 1 

£ £ VijSij{x,t), ( 33 ) 

i=-M m j=-Mt m t = Mt -f Nt + 1. 

The m s • m t unknown coefficients {u^} are determined by orthogonalizing the residual with 
respect to the set of sine functions {S p S;}IZZmV,":X- This y ields the discrete Gherkin 
system 

(£u mmmt -f,s p s;) = o ( 3 -4) 

for p = -M X) • • • t N x and q = -M t , The inner product (•, •) is taken to be 

(F,G) = r f F{x, t)G(x, t)w(x, t)dxdt (3.5) 

Jo J o 

with the weight 

w(x,i) = w(x)w'(i) = W(x)r'H-i(t))-i. (3.6) 
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The expressions (3.1), (3.4), and (3.5) are combined to form the system 

+ Jo Jo S p { x )S*(t)w(x)w*(t)dxdt (3.7) 

= jf £ /(*,<)«#) w(x)w*(t)dxdt 
for p = —M x , ,N X and q = —M t , ■ • • , N t . 

In anticipation of the parameter identification problem which motivates this analysis, the 
term EI{x) in (3.7) is expanded as a linear combination of weighted sine functions with four 
Hermite-like algebraic terms added to accommodate the potentially nonzero function and 
derivative values of El at x = 0 and x = 1. Specifically, this parameter basis is taken to be 
{^k}k=-M. 

r 

b -M m (x), k = —M x 

&-M.+ i(x), k = -M x + 1 

M x ) = | v E (x)S k (x), k = —M x + 2, • • • , N x — 2 (3-8) 

^N m — i(x), k — N x 1 

kw*(x), k — N x . 

Here 5jt(x) = S(k, h x ) o <f>(x ) and the basis weight function vg is 

ve(x ) = tu(x) = [x(l — x)] 5 . (3-9) 

The algebraic boundary basis functions are given by 

6_m.+i(x) = (1 - x) 2 [2x 4- 1], 

^.-i(x) = x 2 [2(l-x) + l], 

&-M.(x) = x(l — x) 2 , 

and 

1>N m (x) = “(I -*)X 2 . 
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(3.10) 


The finite dimensional approximation of El then takes the form 

Va 

EI ma (x)= £ c fc ^ fc (x) 

k=-M m 

where the choices Me = M x and Ne = N x are made so as to guarantee a square spatial 
coefficient matrix. In the forward problem, the coefficients are known whereas in 

the corresponding parameter recovery problem, they are unknown and are determined via 
methods to be discussed in Section 4. 

A quick note should be made concerning the choice of parameter basis and the manner 
of expanding EI ma . The two derivative-interpolating boundary basis functions are added 
so that this expansion of EI mB is the same as that used with cantilever or free boundary 
conditions. The choice of (3.9) for basis weight is certainly sufficient and proves to be 
beneficial when incorporating this forward scheme into a numerical method for solving the 
parameter recovery problem as described in Section 4. 

The expansion (3.10) is substituted into (3.7), and integration by parts is used to transfer 
the derivatives onto the product S p wS*w*. As detailed in [18], the weight choice (3.6) guar- 
antees that all boundary terms vanish. The resulting integrals are evaluated via Theorem 2.2, 
or when possible, Theorem 2.1. The requirement 

\ei(x)u(x,t)\ < Kx a +* (i - xf+f r + h - s , 


where the “homogeneous” part of El is 


£X{x) = EI(x)-EI(0)b_ Mm+1 (x)-EI(l)b„ m _ 1 (x)-Er(0)b_ Mm (x)-Er(l)b Nm (x), (3.11) 

guarantees the decay needed to truncate the infinite quadrature rule as specified by (2.10). 
With a,/3 , 7 and 8 specified and M x chosen, the choices 


h x = 


I ird 
aM’ 


N x = 


h t = h xt 

a 


M x + 1 , 


a 


M e = -M x + 1 , 
7 


( 3.12) 

(3.13) 

(3.14) 

(3.15) 
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and 


N, = []m, + l] (3.16) 

for the stepsizes and summation limits balance the asymptotic errors to at least order 
0 . This rate results from the presence of a sine function in the integral. In the 

above expressions, [•] denotes the greatest integer function. Note that the +1 is unnecessary 
when |M X , ^M x , or jM t are integers. 

In many time-dependent problems, the solution decays exponentially at infinity; that is, 
the solution satisfies 

\€l(x)u(x,t)\ < Kx a +1{ 1 - x) 0+ ir+h- St . (3.17) 

With this supposition, Lund [11] shows that the condition (3.16) can be replaced by 

N t = (jMtht'j + l] . (3.18) 

The selection Nt in (3.18) significantly reduces the size of the discrete system with no loss 
of accuracy. 

Given M x , N x , M t , N t and h = h s = h t as defined above, the discrete system for (3.1) is 

A(EI)UCj + C.UAj = G. (3.19) 


Here 


p 

III 

$* 

(3.20) 

III 

(3.21) 

and 



(3.22) 

where T>(t]) denotes the diagonal matrix with entries r)(x_M m ),' • 

' tV( x N m )- The m x X m ( 


matrices U and F are defined componentwise by 


[^]*j — u »i 


and 

= /(*»>^i)- 
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It should be noted that the ordering of the coefficients in U mimics that used in most 
standard time- differencing schemes. This is a matter of convenience since the Sinc-Galerkin 
method is not bound by any specific ordering of the grid. 

To simplity notation when specifying the spatial and temporal matrices A(EI) and A t 
respectively, let I^\ l = 0,1, 2, 3, 4 denote the matrices whose pt-th entry is Sjp from 
(2.12) - (2.16). As shown in [13], the m t x m t matrix A t is given by 

A t = P((T)*). (3.23) 

The m x x m x matrix A(EI) has the form 

A(EI) = [* (a) 2>(ft<.)) + 2$< 3 >2>(p*(o) + $ (4) £>(P*< o))] . (3.24) 

The notation 7}(p#(o)» ^ = 0,1,2 denotes the diagonal matrices containing the components 
of the vectors 

p* (<) = tf (/) c, * = 0,1,2 (3.25) 

where c= ■ • • , c iv.] T - The matrices j = 2,3,4 and l — 0,1,2 are defined 

componentwise by 

[*% = (3.26) 

and 

[* W U = #(*()• (3.27) 

with the notation on the right-hand sides of (3.26) and (3.27) indicating the j- th and *-th 
derivatives, respectively. 

To illustrate the dependence of j = 2,3,4 and 'i r ^, * = 0,1,2 on fundamental 

matrices, the respective expansions are listed below. The diagonal matrices T) and the 
matrices /W, l = 0,1, 2, 3, 4 have sizes consistent with the following range of indices i,p, 
and k {-M x < t < N x , —M x < p < N Xi -M x + 2 < k < N x - 2). From (3.26) it follows that 
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and 


$ (s > = lzlWV(w<fi') + ^-lWV(2w') + I^V(w") t 
h l h * 

$( 3 ) = -LlWt>(w{<f>’) 2 ) + ^I^V{Zw'(f>' + Zwf) 

+h ,mv ( 3w " +3w ^ +w 9) +,mv {T) 

$( 4 ) = J_/( 4 )-p ( 'w(4 >') 3 ) 4- ^4i u'(<f>') 2 4- 6w<f>'<f >' ') 

+ _L/(2)p ^ 6 u>V + 1 2w'<f>" + 4u/<T 4 3 

+^/Wp + 6»"^ + 4»'^- + w^j 

+/ (0) p ^-). 

V *? / 

For £ = 0, 1,2 the m x x m x matrices in (3.27) are given by 


$(0 = 


m : m : nil) : r(0 : £(') 

® -M B • °-A<r.+l • ^ • “N.-l • °N m 


(3.28) 


(3.29) 


(3.30) 


(3.31) 


where b]P = [6j^(®_*f.), • ■ • , ^k\ x N m )] T for k = —M x , —M x + 1 , N x — 1 and N x . Again, the 
superscript l indicates the ^-th derivative. The m x X (m x — 4) matrices are 


= V(v E )l(°\ (3.32) 

BW = ~ V{v E <f>')I (1) + V{v E )lW (3.33) 

and 

B m = * V [veWY) /<» - + 2»W')/« + X>K)/ (0) . (3.34) 

The negative signs that appear in the definitions of B ^ and B ^ result from the transposing 
of jW. Again, it is noted that in (3.32) - (3.34), the m x x (m x — 4) matrices I^\l = 0,1,2 
have components ^ as defined in (2.12) - (2.14). 

Various methods exist for solving matrix systems of the form (3.19), one of which derives 
from the generalized Schur decomposition (page 396 of [4]). As guaranteed by the results of 
Moler and Stewart [15], there exist unitary matrices Qi,Zi,Q 2 and Z 2 such that 
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Q[A x Zi = P 
Q[C X Z i = R 
Q* 2 C t Z 2 = S 


Q 2 AtZ 2 = T 

where P } R,S, and T are upper triangular. If Y = ZIUZ 2 and C = Q[GQ 2 , then (3.19) 
transforms to 


PYT* + RYS' = C. 


By comparing the A:-th columns, one finds that 

n n 

p Y + R - Ck 

i=k j=k 

which yields 

n n 

(t kk P + s kk R)y k = c fc - P Y, t ki yj - R Y Sk jVi (3.35) 

j=fc+l i=k + 1 

(for convenience, it is assumed that all matrices are n X n and indexed from 1 to n). With 
the assumption that the matrix ( t kk P + SkkR) is nonsingular, the solution to (3.35) is easily 
found by recursively solving triangular systems. 

Although this algorithm does require complex algebra, it is both robust and efficient 
and requires no assumptions concerning the diagonalizability of the component matrices. 
It should be noted that a “real” version of this algorithm also exists [3]. In this latter 
algorithm, Q\, Z\,Q 2 , and Z 2 are orthogonal with P, S quasi-upper triangular and R f T 
upper triangular. 

A Sinc-Galerkin method for the more general problem (3.2) can be developed in a similar 
manner once a suitable basis has been determined for discretizing the spatial variable. To 
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this end, define the set of spatial basis function by 


CM = 


B_ Mm _ A (x), 

•xf 

1 

sf 

1 

II 

B_ Mm _ 3 (x), 

i = -M x - 3 

B- Af.-j(x), 

i = —M a — 2 

B-m.~ i(z), 

♦*. 

II 

1 

* 

1 

b-* 

v(x)Si(x), 

* = -M x , • • • , 

Bn.+i(x), 

i = N x + 1 

B N.+2(x), 

i = N x + 2 

B Ar.+3(z), 

i = N x + 3 

B N.+i(x), 

i = N x + 4. 


Here 5,-(x) = S(i, h x ) o <f>[x ) and the basis weight v(x) is taken to be 

u(x) = [x(l — as)] 3 . 

The boundary basis functions are 

B_ m .-i(x) = (1 - x) 4 [20x 3 + 10x 2 + 4x + 1], 

B Nm+1 (x) = x 4 [20(l - x) 3 + 10(1 - x) 2 + 4(1 - x) + 1], 
B- Af.-j(x) = x(l - x) 4 [10x 2 + 4x -f 1], 
B Nm + 2 (x) = -x 4 (l - x)[10(l - x) 2 + 4(1 - x) + 1], 
B_m.-z{x) = x 2 (l - x) 4 j^2x + i] , 

B„.„(x) = x*(l - x) J [ 2(1 - x) + i] , 


and 


B-m .- <(*) = g® 3 (l ~ x ) 4 » 


B N m +i{x) = - *) 3 . 


(3.36) 


(3.37) 
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A brief note concerning the choice of spatial basis is in order at this point. First, since 
^jj[5(t, & x ) o <^(x)] , l = 1, 2, • • ■ , are undefined at x = 0 and * = 1, Borne basis modifica- 
tions must be made when solving problems with nonzero boundary conditions (see also the 
definition of in (3.8)). With the cantilever boundary conditions of (3.2), it is tempting to 
use fewer algebraic boundary basis functions and the basis weight 

v(x) = x(l — x) 3 


but in many problems this results in nonzero boundary terms when integrating by parts. 
By using the symmetric basis weight v(x) = [x(l — x)] 3 and a full complement of algebraic 
terms, this pitfall can be avoided. Furthermore, the spatial basis {£,} as defined in (3.36) 
can be used for problems with free boundary conditions, thus providing consistency to the 
method. 

The basis for the problem (3.2) is then taken to be {(iSj} with Sj(t) = S(j,h t ) o T(f), 
and the approximate solution is defined to be 


«m.m f (®,0 = 2 X) u <i&( x ) 5 j(0 

i= — j= — Mt 

N t 

+ ]£ S J *(t){U-Af.- 3 ,jC-M 11 -3 (*)+« -M.- 4,jC — 4 (x) 

j=-M t 

+UW. + l,jGv.+l(®) + U|V.+2,jCtf.+2(*)} 


+{a(0C-M.-i(x) +/5(t)C-M.-2(z) + 7( t )Cw.+3(z) + 5(<)C W . + 4(x)} 

(3.38) 


where 


7(0 = 


1 

El( 1) 


7(0 


and 

= EI(1)^ ~ [£/(l)] 2 ^' 

The functions 7(f) and 5(0 are well-defined since EI(x ) is assumed positive on [0,1]. It 
should be noted that the approximate solution does satisfy the boundary conditions in (3.2). 

The (m x + 4) • m t unknowns {«,,} in (3.38) are determined by orthogonalizing the resid- 
ual with respect to the sine functions {<S’ p (x)S'*(0}p=-M,-v , . > Ar B +3* This Petrov-Galerkin 
approach is in contrast to those Galerkin methods in which the residual is orthogonalized 
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with respect to the basis and is done to take advantage of the exponential accuracy of point 
evaluation in the quadrature. This yields the discrete system 


(£u m „ m , — f, S p S*) — 0, 


-M t — 2<p<N x + 2 


—M t <q< N t 


where (■, •) is defined in (3.5) with w(x,t) = (T(<))“i. 

Appropriate integration by parts and application of the sine quadrature rules, as discussed 
in [18], yields the matrix equation 


A(E1)UC? + C x MUAj = G 

where Ct,C x , and At are defined in (3.20), (3.21), and (3.23) respectively, and 

The ( m x + 4) X m x matrices U and are defined componentwise by 

Wii = 

and 

F]ij = 7(*i, tj) 

where 

70M) = /(*,<) - 

-CWt)B Na+3 (x)) - C(6(t)B Nm+4 (x)). 

The (m c + 4) x ( m x + 4) matrix M has the form 


(3.39) 



I 1 

1 |0 

l 

1 

M = 

. 1.1 
1>L2 1 bn | Dm | 

, 1 . 

1 Oftl 1 b R2 


1 1 1 

1 l 0 1 

1 

1 


where the m x x m x submatrix Dm and the (m, +4) x 1 vectors are given by 

Dm = D(v), 
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1>L2 = 4 ) 1 , 

hi = 3)1, 

hi = V(Bif m+ i)l 

and 

h 2 = E>(Bff m+2 )l. 

Here 1 is simply an (m s + 4) X 1 vector of ones. 

The spatial matrix A(EI ) can be constructed as follows. Let and p${i) be defined 

as they were in (3.26), (3.27), and (3.25), respectively (with / = 0,1,2 and j = 2,3,4). Note 
that in the definitions now, the index ranges are — M x <i< N x , —M x — 2<p<N x + 2, and 
—M x —2 < k < N x + 2, and the (m x + 4)x 1 coefficient vector is now c = [c_m.-2, • • • » c„.+j] . 
Hence 'pW and p* ( o have the sizes (m, + 4) x m x , m x x (m x + 4) and m x x 1, respectively. 
Furthermore, let denote the (Tn x + 4) x m x matrix which is defined componentwise by 

and let c = [ c_m ■ ■ • , cjv.] t . Finally, fori = — M x — 4, • • • , —M x — 1 andi = N x + 1, • ■ - ,N X + 4, 
let Oj denote the (m x + 4) x 1 vectors 

Bi = $"v{v e b';)c + v ^(£; c fl,")"j 1] . 

Here 

EI e (x ) = c_Af._i&_M.-i(*) + c_M.- 2 &-M.-a(a:) + c//. + i&Ar. + i(*) 4- CAr„ +a &jv. + 3(®) 

and 1 is simply the m x x 1 vector of ones. The (m x + 4) x (m x + 4) matrix A(EI) is then 

A(EI) = [a_M m ~ 4 : B-M.-3 ' A m : ajv.+i : ^.+ 2 ] (3.40) 

where the (m x + 4) x m x submatrix A m is given by 

= [$ (2) £>(v)P(p*<>>) + 2$Wv(v)V(p* {l) ) + $< 4 >£>(v)27(p* ( o))] . 

It should be noted that the coefficient matrix A(EI) in (3.40) differs from that arising 
in the fixed boundary problem, (3.24), only in the presence of the diagonal multipliers V(v) 
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and the addition of border vectors. Hence the method is easily adapted when changing the 
boundary conditions. Furthermore, the matrices 4", and 4^ can be expanded in terms 
of fundamental matrices in a manner similar to that in (3.28) - (3.34), thus simplifying the 
implementation of the method. 

With A(EI), A t , C x , C t , M, and G thus specified, the system (3.39) can be solved via the 
generalized Schur algorithm (3.35). 

The final implementation issue for the forward problem (3.2) concerns the choice of decay 
parameters a, (3, 7, and 6. As discussed in [18], the weight choice w(x,t) = (T(t)) _ * yields 
the decay condition 

\£l{x)U{x,t)\ < x“+ 3 (l - + 3 r+ 5 e -“ (3.41) 

where £I(x) is defined in (3.11) and Z/(x,t) is that part of the true solution which is ap- 
proximated by 

«a(M) = £ ViiCiWSji*) 

i=~M m j=~ Mt 

(this can be formally obtained by subtracting all boundary contributions from the true 
solution u(x, t)). With the decay parameters specified and M m chosen, the remaining stepsizes 
and summation limits are given by (3.12) - (3.16). 


4 The Finite Dimensional Minimization Problem 

As noted in the introduction, the minimization problem 

min TJEI) 

EitQ v ' 

where 

X.(EI) = ^{\\IC(EI)-df+ a \\Em1 1 } 

is infinite dimensional and thus must be replaced by a sequence of finite dimensional problems 
before a viable numerical scheme can be developed. Following from (3.10), the approximating 
admissible parameter sets are taken to be 

Q„, = [ El„, ■. EI„ a (x) = ^ c ,iMx) 

I 


| j m E ~ Me + Ne + 1 
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with the basis {V’fc} defined in (3.8). The summation limits depend on the boundary condi- 
tions with Me = M x , Ng = N x for fixed boundary conditions and Me = M x - f 2, Ne = N x + 2 
for cantilever boundary conditions. The associated finite dimensional optimization problem 
can then be formulated as 

min f.(EI m .) (4.1) 

where 

= i{||A-(E7 m ,) - d'll 1 + or||B/„,||J}. (4.2) 

Note that in solving the minimization problem (4.1), one is actually solving for the vector 

_ A 

c = [c-M a r • • , cn b ] € ET 1 " which minimizes T a . 

With n p and n, specifying the number of spatial and temporal observation points, re- 
spectively, the approximation : 2R m * — ► ]R np ' n * to fC(EI) is obtained by applying 

a discrete analogue of the observation operator C in (1.6) to Um amt in (3.3) or (3.38). If the 
set of observation points {(x p ,i,)} p ^}j"l£* can represented as a tensor product of spatial 
and temporal points, then K(EI m3 ) has the representation 

K(EI mB ) = C co{U) (4.3) 

where the matrix U solves either (3.19) or (3.39) depending on the boundary conditions. For 
the fixed boundary problem (3.1), the matrix concatenation co{U ) is the vector in JRT' m ' mi 
which is obtained by successively stacking the columns of the m x x m t matrix U. C is an 
(n p -rig) x (m, -m t ) evaluation matrix which can be formulated as follows. Define the n p x m x 
spatial evaluation matrix E x to have components 


[E x ] Pti = Si{x p ), 1 <p<n p , -M x <i<N x 
and let the n q x m t temporal evaluation matrix E t have the components 


Then 


C = E t ® E x . 


Similar formulations for C and co(C7) can be used if the cantilever boundary value problem 
is being considered. It is noted that if the set of observation points is not rectangular as 
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described above, then point evaluation can be done directly via (3.3) or (3.38). This latter 
option is less efficient, however, than that defined in (4.3). 

The discrete penalty functional ||2?/ ma ||g is formed by substituting the expansion (3.10) 
into the definition (1.8). This yields 

l|£A-X = + £ » l T QS 

JO Jo 

where the m E x m E matrix Q = Qd Q / has components 

[Qd]kt « ( -M e <k,i<N E 

Jo 

and 

[Qf}ki ~ / rpk{x)rpt{x)dx, —Me <k,l< N E . 

Jo 

The matrix entries are approximations in the sense that sine quadrature rules are used to 
evaluate many of the integrals. 

For the choice of basis functions in (3.8), the matrix Qd is given by 

4 6 0 -6 2 

6 12 0 -12 6 

-» -* « -* -* 

0 0 Qd 0 0 

-6 -12 0 12 -6 

26 0 -64 

Integration by parts and the application of the sine quadrature formula (2.11) yields the 
(m E — 4) x (mjs — 4) submatrix 




«'-E iW -!E jW + IS jW 


where again, 1^,1 = 0,2,4, denote the matrices whose pt-th entries are 6^ from (2.12), 
(2.14), and (2.16), respectively. The zeroing of all other quadrature terms is a result of the 
choice ve(x) = (<f>'(x))~l = [x(l — x)]$. The notation 0 simply denotes a zero vector of 
length m E — 4. Because 1^ is positive definite and is negative definite (see [17]), the 
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matrix Qd is nonnegative definite with zero eigenvalues resulting from the configuration of 
the outer four columns. 

Direct integration and sine quadrature are used to obtain the matrix 

I 11 -T 13 -1 

105 310 " t2 420 140 

II 13 -T JL _LL 

210 35 ’ *1 70 420 

qi2 qtl Qf 9rl 9r2 

13 _9_ -*T 13 11 

420 70 " rl 35 120 

-1 13 r*T 11 1 

140 420 " r2 120 105 

Here 

Q f = h x V{x 3 { 1-x) 3 ) 

where denotes the diagonal matrix with entries i/(®-Af B +2)» * • • > v{ x Nb- 2 )- R eca ^ that 

the sine gridpoints are defined in (2.7). The vectors qti,qi 2 ,q r i and q r2 have components 

[qn]k = h x x k (l - *fc) 3 [2*fc + 1], 

[<f«]k = h x (l - ® fc )ifc[2(l — ifc) + 1], 

[9rl]fe = h x xl(l — Xfc) 3 , 

and 

[5-2]fc = -h x (l - Xk) 2 xl 

for k = —Me + 2, • • • , Ne — 2. The matrix Qf is strictly positive definite. 

Although the matrix Q is full, it is very efficient to construct since the Toeplitz matrices 
jf°), jW and I W are also needed in the forward solver. For e > 0, Q is symmetric and positive 
definite and hence has a Cholesky decomposition Q = R T R where R is upper triangular. It 
then follows that the penalty term ||i?J ma ||Q yields the quadratic form 

c t R t Rc= ||ite|| a (4.4) 

where || • || denotes the Euclidean norm. This representation for the penalty functional is 
particularly useful both when implementing a scheme to solve the minimization problem and 
when plotting the .L-curve to determine a suitable regularization parameter oc (see [6]). 
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To highlight the dependence of the functional T a in (4.2) on the unknown vector 
C = [c_ Mg » ‘ ' i ) let 


K{c) = k{EI mm ) = C co{U (c?)) 

where U{c ) solves either (3.19) or (3.39) for a given expansion El ma (x) = 
Noting (4.4), the optimization problem (4.1) can be replaced by 


X) c ki>k(x). 
k=-M a 


min 


r.(*) 


(4-5) 


where 

T.W = i{||itr(?) - iff + »||KTf }• (4.6) 

To obtain a minimizer for the nonlinear functional T a (c), a quasi-Newton/trust region iter- 
ation [2] 

c k+l — c k s k 


is used. Here s* solves the quadratic programming problem 


,™Sl B \ + K'fa)*- rf l Z + “11^(4 + a )|| 2 } 

subject to ||s || < Sk with K'(ck) denoting the Jacobian of K at c* . The trust region radius 
8k is chosen so that T a (c ) has sufficient decrease at each iteration to guarantee convergence 
to a local minimizer of T a (for further details about the theory and implementation of the 
trust region algorithm, see [2] or [8]). An important numerical issue in the implementation 
of the trust region scheme is the formulation of the derivative of the operator K. Here the 
derivative, or Jacobian, is an ( n p • n 9 ) x tre matrix whose v-th column is given by 


[*'(*)]. = Jim j[K{c + Te„) - *(<?)] 

where the standard unit vector e„ has components 


[K]k = 8„k = < 


1 

0 


if As = t', — Me < k Ne 

otherwise. 


In the examples of the next section, the Jacobians were calculated with a standard forward 
difference scheme. This scheme is easy to implement and accurate enough for the purposes 
of the method. If further efficiency is desired, a directional derivative scheme such as that 
described in [12] can be used. 
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5 Implementation and Numerical Examples 


The four examples reported in this section were selected from a large collection of problems 
to which the Sinc-Galerkin method was applied. The results are representative of those 
obtained for other problems. 

The first example demonstrates the application of the Sinc-Galerkin method to a model 
problem with fixed boundary conditions in which the state solution was sampled directly; 
that is, no external noise was added to the data. To demonstrate the feasibility of the method 
for problems with noisy data, the same problem is revisited in Example 5.2 but with pseudo- 
random white noise added to the data. In the third example, the parameter to be recovered 
is the shifted Gaussian function that was discussed for second-order examples in [12] and [19]. 
Again, pseudo-random white noise is added to the data. The final example demonstrates 
the application of the method to a problem with cantilever boundary conditions as modeled 
by (1.2). Hence in this example, the parameter El appears both in the spatial operator and 
in the boundary conditions. 

In all examples, d = f (see (2.4) and (2.6)). The errors for the method are reported on 
the set of uniform gridpoints 

U = { Z 0y Zi , • ■ • ,Zioo} . 

(5.1) 

Zk = kt ) k = 0, 1, * ■ • , 100 

with stepsize i = ^ With El and EI ma denoting the true and approximate parameters 
respectively, the errors are reported as 

IIMOII = I (52) 

The error results are tabulated in the form .ana - 7 which represents .ana x 10 - ' r and all 
problems were run with sixteen place accuracy on a Vax 8550. 

A very important practical consideration is the choice of the regularization parameter a 
for a given (error-contaminated) data set. If the error in the data is discrete and random, 
then under certain conditions the method of Generalized Cross Validation (GCV) can be used 
to determine a suitable value of a [25]. A second method for determining the regularization 
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parameter is to plot the norm of the penalty functional, ||i?c a ||, versus the norm of the 
residual, [|/ i f(c 0 ) — d || (see [6] or [14]). Here c a denotes the solution to (4.5). In this way, 
one can qualitatively get an idea of the compromise between the minimization of these two 
quantities. The scheme for determining the “optimal” regularization parameter consists of 
finding those values of a such that (||.ff(c a ) — d ||, ||i2c a ||) lies in the “corner” of the resulting 
curve, known as the Z-curve. The use of this technique for determining suitable choices for 
the regularization parameter is demonstrated in the examples. 

In all four examples, the data was sampled on a regular grid {(x p , <,)} C (0, 1) x (0,2]. 
Nineteen equally spaced points x p = pAx, Ax = .05, were taken in space and four equally 
spaced temporal points t q = qAt, At = .5, were taken for a total of n = 76 data points. In 
all examples, the x 1 initial vector cq = [.5, .5, .5, • • • , .5, .5, .5] r was used. 

Finally, it should be noted that in the examples, the symbol a is used to denote both 
the regularization parameter and the sine decay parameter. The use of this symbol for both 
quantities is well established in the literature and is thus difficult to avoid in this setting. 
It should be obvious from the context however, which quantity is being discussed and there 
should be no ambiguity resulting from the dual use of this symbol. 

Example 5.1. 

d 7 u d 2 / c? 2 u\ 

W + d^[ EHx) dJ j = /(x ' () ' (>0 0<I<1 

tt(0, t ) = u(l, t) = 0, t > 0 

^(o,0=|(M) = o, <>o 

Qxi 

u(x, 0) = 0) = 0, 0 < X < 1 

The forcing function f(x, t ) is consistent with the true stiffness parameter EI(x ) = l+sin( 7 rx) 
and the state solution u(x, i) = x(l — ®) sin(47rx)t 2 e -< . For these functions, the choices 
a = P = 7=f and 5=1 satisfy the decay condition (3.17). No noise was added so 
the data consisted of direct measurements of the state solution. For varying values of the 
regularization parameter a, the L-curve is plotted in Figure 3. Note that the value a = 10 -8 
yields a point (||.fiT(c a ) — <f||, ||.ffc a ||) in the “corner” of the curve. The uniform errors for 
a = 10 -8 are reported in Table 1 with the first four columns indicating the index limits 
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for the expansion of the state variable and the fifth column indicating the number of basis 
functions used in the expansion of EI mB . The convergence of the method is demonstrated 
both by the results in the last column of Table 1 and by Figure 4 which shows the true and 
approximate stiffness parameters with a = 10 -8 . 


M x 

N x 

M, 

N t 

m E 

\\Eium 

8 

8 

8 

4 

17 

.1869 - 0 

16 

16 

16 

6 

33 

.2482 - 1 

24 

24 

24 

7 

49 

.1463 - 2 



1.5 2 2.5 3 3.5 4 

Reaidu&l Norm — «f || xlO' 1 

Figure 3. The Tikhonov L-Curve for Example 5.1. 
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Figure 4. True and Approximate Stiffness Parameters for Example 5.1 with a = 10 8 
- • “ = 17), (m E = 33), • • • (m E = 49), (True). 


Example 5.2. 

In this example, the true parameter and state solution are the same as those in Example 5.1, 

and hence EI(x) = 1 -f sin(Trx) and u(x,t ) = x(l — x) sin(47rx)f J e -( . To the data however, 

we added a pseudo-random noise vector e from a Gaussian distribution with mean 0 and 

— * 

standard deviation a chosen so that the noise- to-signal ratio a/||d || = .01 ; that is, noise 
= 1% of the signal. The L-curve is plotted in Figure 5. Note that the values a = 10~ 5 and 
a = 5 x 10~ 6 yield points (||/f(o,) — c?||, ||7Zc a ||) in the “corner” of the curve. For m E = 33, 
the uniform errors obtained with a = 10 -3 , a = 10 -5 and a = 10 -1 ° are given in Table 2. 
Corresponding plots of the true and approximate stiffness parameters are shown in Figure 6. 
Note that the “corner” value a = 10 -5 provides a very good choice for the regularization 
parameter whereas a = 10 -10 is not large enough to damp out the error contributions due 
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to the smaller singular values. Finally, the choice a = 10 -3 causes too much smoothing and 
information about the parameter is lost. The results from this example demonstrate the 
viability of the method for problems with noisy data. 



a = 10" 3 a = 10" 5 a = 10~ 10 


.8503 - 0 .2228 - 1 .3840 - 1 


Table 2. Errors on the Uniform Grid U with mg = 33 in Example 5.1 



Figure 5. The Tikhonov L-Curve for Example 5.2. 






x-axis 


Figure 6. True and Approximate Stiffness Parameters for Example 5.2 with m r = 33 
(a = 10- 3 ), (a = 10- 5 ), •■•(« = 10- 10 ), (True). 

Example 5.8. 


i>0 0<I<1 
u(0, t) = v(l, <) = 0, i> 0 

§>0=^(M) = o, <>o 

u(x, 0) = ^(x,0) = °> 0 < X < 1 

In this example the stiffness parameter to be recovered is the shifted Gaussian function 
EI(x ) = 1 + The state solution u(x,t) = sin 3 (xx)f a e~* yields the decay param- 

eters a = (3 = 7 = f and 8 = 1 as dictated by (3.17). Pseudo-random noise is again added 
to the data in the manner described in the last example. As seen in Figure 7, the Tikhonov 
parameter values a = 10~ 6 through a — 5 x 10 8 yield points (||A(ca) — <T||, ||i2c a ||) 
in the “corner” of the L-curve. For the “corner” value a = 10" 7 , numerical results with 
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mg = 17, mg; = 33 and mg = 49 are reported in Table 3 and Figure 8. In spite of the noise 
in the data, both sets of results demonstrate that the method converges as the number of 
basis functions is increased. 


M x 

N x 

M, 

N t 

m E 

IIEMOII 

8 

8 

8 

4 

17 

.1811-0 

16 

16 

16 

6 

33 

.4191 - 1 

24 

24 

24 

7 

49 

.1058 - 1 


Table 3. Errors on the Uniform Grid U with a = 10 7 in Example 5.3. 



Figure 7. The Tikhonov L-Curve for Example 5.3. 
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Figure 8. True and Approximate Stiffness Parameters for Example 5.3 with a = 10 7 
- • - (m E = 33), (m E = 49), (True). 

Example 5-4- 


d 2 u d 2 d 2 u\ 

ai 7 + fe 7 j = ^ x,t ^ t>0 0<I<1 

u(0,i) = 0, fE/~J(l,t) = 0, t> 0 
u(x,0) = ^(x,0) = 0, 0 < x < 1. 

This example demonstrates the Sinc-Galerkin method for a problem with cantilever boundary 
conditions and hence the stiffness parameter appears both in the spatial operator and in the 
boundary conditions themselves. The true stiffness parameter is EI(x) = l+sin(7rx) and the 
state solution is u(x,t) — sin 3 (7rx)f 2 e -< . For these functions, the choices a = /? = l, 7 =y 
and 6=1 satisfy the decay condition (3.41). No noise was added so the data consisted 
of direct measurements of the state solution. Since the L-curve was very similar to that of 
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Example 5.1, the regularization parameter was taken to be a = 10~ 8 . The uniform errors 
for this choice are reported in Table 4 and the true and approximate parameters are shown 
in Figure 9. When comparing Tables 4 and 1, it is noted that the index limits differ as a 
result of the choice ot — (3 = 1 for the decay parameters. The smaller values of ot and /3 also 
indicate why the errors here are slightly larger than those in Example 5.1. Finally, because 
M e = M x + 2 and N E — N x + 2 for these boundary conditions, the number of basis functions, 
m E} used in the expansion of EI rrXB also differs from the number used in Example 5.1 where 
M e = M x and N E = N x . Both the table and the figure demonstrate the convergence of the 
method for problems with cantilever boundary conditions. 


M x 

N x 

M t 

N t 

m E 


8 

8 

6 

3 

21 

.1589 - 0 

16 

16 

11 

4 

37 

.2589 - 1 

24 

24 

16 

6 

53 

.1334 -1 


Table 4. Errors on the Uniform Grid U with a = 10 8 in Example 5.4. 



Figure 9. True and Approximate Stiffness Parameters for Example 5.4 with a - 10 8 
(m E = 21), (m B = 37), ■ • ■ (m E = 53), (True). 
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